Using basis sets of scar functions 
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We present a method to efficiently compute the eigenfunctions of classically chaotic systems. The 
key point is the definition of a modified Gram-Schmidt procedure which selects the most suitable 
elements from a basis set of scar functions localized along the shortest periodic orbits of the system. 
In this way, one benefits from the semiclassical dynamical properties of such functions. The perfor- 
mance of the method is assessed by presenting an application to a quartic two dimensional oscillator 
whose classical dynamics are highly chaotic. We have been able to compute the eigenfunctions of the 
system using a small basis set. An estimate of the basis size is obtained from the mean participation 
ratio. A thorough analysis of the results using different indicators, such as eigenstate reconstruction 
in the local representation, scar intensities, participation ratios, and error bounds, is also presented. 



I. INTRODUCTION 

The vast majority of methods to obtain quantum sta- 
tionary states rely on the expansion of the correspond- 
ing wave functions in a basis set of suitable basis func- 
tions that can be made (approximately) complete, on 
which the Hamiltonian of the system is diagonalized. The 
choice of the basis set is then critical for the efficiency of 
the method. This issue is particularly important in the 
case of heavy particle dynamics or in the semiclassical 
limit, where these functions oscillate considerably. The 
situation is even worse for very chaotic or ergodic sys- 
tems, as those in which we are interested in this paper. 

Several methods have been proposed in the literature. 
The simplest procedure uses products of harmonic os- 
cillator eigenfunctions, something which works well to 
describe a good number of low- lying states, but gets pro- 
gressively poor as energy increases due to anharmonici- 
ties (see, for example, Refs. [T]and[2]). Going to the other 
extreme, other methods have been proposed making use 
of the semiclassical information derived from quantized 
invariant classical structures [3] , that render excellent re- 
sults HE]. 

In this paper we investigate the feasibility of using scar 
functions, localized over short periodic orbits (POs), as 
a basis set to efficiently compute the eigenstates of clas- 
sically chaotic Hamiltonian systems. 

The term "scar" was introduced by Heller in a sem- 
inal paper [6 to describe the dramatic enhancement of 
quantum probability density that takes place along POs 
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in some eigenfunctions of the Bunimovich stadium bil- 
liard, as a result of the recurrences along the scarring 
orbit. The relevance of unstable POs in the quantiza- 
tion of classically chaotic systems had been previously 
pointed out by Gutzwiller in his celebrated trace formula 
(GTF) [7j- Other fundamental contributions to the the- 
ory of scars [8] were made by Bogomolny [9 , who showed 
how this extra density is obtained by averaging in config- 
uration space groups of eigenfunctions in an energy win- 
dow around Bohr-Sommerfeld (BS) quantized energies in 
the h —> limit. The corresponding phase space version 
using Wigner functions was investigated by Berry [TO] . 
Other interesting aspects of scarring, such as the role of 
homoclinic and heteroclinic quantized circuits [TTJ [12], 
the influence of bifurcations (in systems with mixed dy- 
namics) [13], the scarring of individual resonance eigen- 
states in open systems [14], or relativistic scarring [15] 
have also been discussed in the literature. Scars have also 
been experimentally observed in many different contexts, 
including microwave cavities [16] , semiconductor nanode- 
vices [17 , optical microcavitities [18], optical fibers [19], 
and graphene sheets [20] . 

Different methods have been described in the literature 
to systematically construct functions localized on unsta- 
ble POs (hereafter called scar functions). Polavieja et 
al. averaged groups of eigenstates using the short-time 
true quantum dynamics of the system [2T] , Vergini and 
coworkers used the short POs theory [22] and obtained 
scar functions by combination of resonances of POs over 
which condition of minimum energy dispersion is im- 
posed, thus including the semiclassical dynamics around 
the scarring PO up to the Ehrenfest time [23] . Sibert 
et al. [24 and Revuelta et al. [25] extended the method 
to smooth potential systems. Also, Vagog et al. extend 
to unstable POs the asymptotic boundary layer method 
to calculate stable microresonator localized modes [26] . 
These scar functions appear not only well localized in 
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configuration and phase space, but they also present a 
very low dispersion in energy [27], and this property 
makes of them good candidates a priori to form an ef- 
ficient basis set for the calculation of the eigenstates of 
classically chaotic systems. An additional advantage of 
using this kind of basis functions, which are based on dy- 
namical information, is that they should allow an easy 
and straightforward identification of the underlying in- 
variant classical structures that are relevant for the semi- 
classical description of individual states of a chaotic sys- 
tem. 

In this paper we introduce a new method to construct 
basis sets formed by the scar functions described be- 
fore [23-25 that can be used to efficiently compute the 
eigenvalues and eigenfunctions of classically chaotic sys- 
tems with smooth potentials. 

This methods exploits a simple semiclassical idea, 
based on the well known Weyl law for closed systems, 
which gives an intuitive explanation of how the quantum 
states of a system "fill" the corresponding phase space [3]. 
Put in numerical terms, the associated volume divided by 
that of a Planck cell (that taken by a single state) gives 
a semiclassical estimation of the generated Hilbert space 
size, 
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where d is the dimensionality of the problem, and O the 
Heaviside function. The application of this prescription, 
i.e. the calculation of the phase space integral in the 
above expression, is particularly simple when one aims 
at calculating all states up to a given energy, E , which 
then provides an easy way to compute a minimum bound 
to the dimension of the required basis set. Obviously, it 
is always advisable to increase this number a little bit 
to account for the border effects, in order to obtain a 
better description of the states localized on this region 
of phase space. This type of strategy has long been pro- 
posed in the literature. For example, Heller et al. [5] 
choose to fill the relevant phase space up to the consid- 
ered energy with coherent states or linear superpositions 
of such states placed along quantized trajectories. No- 
tice that the usual basis sets, constructed for example 
with (orthogonal) products of harmonic oscillator func- 
tions on each coordinate, are less efficient since they are 
worse adapted to the relevant phase space, unnecessarily 
extending into the classically forbidden regions. Bogo- 
molny [4 used semiclassical functions distributed in a 
narrow crust around a given energy shell, thus covering 
a phase space volume given by the surface of the mean 
energy shell times the energy width of the functions. Our 
method is similar in spirit to this one, and the required 
phase space up to a given energy is then filled up with 
a succession of overlying layers, one in top of each other, 
in an onion-like fashion. Very recently [28] , basis sets 
of scar functions defined over short POs have been used 
to compute the eigenstates of the evolution operator in 
open quantum maps. In addition to showing a great per- 



formance for this task, they have proven a very powerful 
tool for the analysis of the behavior of these kind of sys- 
tems. The fact that our scar functions are defined with 
a very low energy dispersion makes that they fill very 
effectively, i.e. with smaller basis sizes, the phase space. 
Notice that some complications arise when constructing 
the basis elements, due to the overlap existing among the 
scar functions. Finally, let us remark that the kind of ba- 
sis sets proposed here can be considered as optimal from 
a semiclassical point of view, since they minimize the dis- 
persion by making a time evolution until Ehrenfest time 
[see Eq. ([l2]) below]. 



The performance of the method is illustrated with an 
application to a highly chaotic coupled quartic oscillator 
with two degrees of freedom that has been extensively 
studied in connection with quantum chaos [TJ [2j [29] . 
We show how our method is able to accurately compute 
the ~2400 low-lying eigenfunctions of the quartic oscil- 
lator using a basis set consisting only of ~2500 elements 
constructed over 18 POs of the system. Furthermore, 
we demonstrate that the method can be used to com- 
pute the eigenfunctons of the system in a small energy 
window using a basis set whose size is of the same or- 
der of magnitude as the number of computed eigenfunc- 
tions. An estimate of this basis size is obtained from 
the mean participation ratio, and it results much smaller 
than those used in other methods. The extension to sys- 
tems of higher dimensionality is straightforward. 



The organization of the paper is as follows. In Sect. [TT] 
we introduce the system that has been chosen to study. 



Section [111] is devoted to the description of the method 
that we have developed, which is based on two central or 
key pillars. First, we use a general procedure able to con- 
struct localized functions along unstable POs for Hamil- 
tonian systems with smooth potentials. Second, we use a 
selective Gram-Schmidt method (SGSM) to select a lin- 
early independent scar function subset from an overcom- 
plete set within a given energy window, thus obtaining 
by direct diagonalization of the corresponding Hamilto- 
nian matrix the desired eigenenergies and eigenfunctions 
with a great degree of accuracy using standard routines. 
We also describe in this section the different mathemat- 
ical tools that will be later used in the analysis of the 
quality of our results. They are: 1) local representation 
functions, 2) scar intensities or contribution of each PO 
to the emergence of an individual eigenfunction, and 3) 
participation ratios, from which an approximated idea of 
the number of basis elements needed to reconstruct an 
eigenfunction can be obtained. In Sect. IV we present 



and discuss the results obtained in the calculation of the 
eigenstates of the system. Finally, in Sect. [V] we sum- 
marize the main conclusions of our work and make some 
final remarks. 
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Figure 1. (color online) (a) Equipotential line for the quartic 
oscillator corresponding to £ = 1. 

(b) y — with P y > Poincare surface of section for a typical 
trajectory of the system at £ = 1. No signs of regular motion 
are apparent. 



II. SYSTEM 

Our system consist of a particle of unit mass moving 
in a quartic potential on the x — y plane 



1-L(P x ,Py,x,y) = 



p2 

y 



^f + ^ 4 + y% (i) 



with the parameter (3 = 1/100. This Hamiltonian has 
been very often used in studies concerning quantum 
chaos [TJ [2j [29j [30]. A plot of the equipotential line 
£ = 1 is shown in Fig. [lja). The corresponding dy- 
namics are highly chaotic; notice that no signs of in- 
variant tori can be identified in Fig. [ljb) , where the 
{y = 0, P y > 0} Poincare surface of section (SOS) for 
a typical trajectory at £ = 1 is shown. However, some 
stable POs exist [3lJ[32], although the area of their sta- 
bility regions are negligible for all practical purposes. 
Another interesting property of Hamiltonian is that, 
due to the fact that the potential is homogeneous, it is 
mechanically similar. This implies that any trajectory, 
(xt,yt, P x ,t, Py,t)i a t a given energy, £, can be scaled to 
another, (x' t ,, y' t ,,P' x t ,,Py t ,), at a different energy, £ by 
using the simple scaling relations 



Vt 
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Vt, Py,t> = 



£' 



1/2 



£ 

£^ 1/2 
£ 



(2) 



Py,t, 



where t' = (£ ' /£ )~ 1 ^H. Notice, that the above relations 
imply that the structure of the phase space is the same 
for all values of the energy. Furthermore, from these ex- 
pressions it can be easily shown that the classical action, 

S t = Jq dr [Px yT + Py,r\-> sca l es as 



3/4 



St- 



PO 



(3) 
fulfills 



Notice that the period of 
T' = (3/4)S T /£, since T = dS/dE. 

Finally, we present in Fig. [2] the eighteen POs for the 
system described by Hamiltonian that will be used 
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Figure 2. Unstable periodic orbits for the quartic oscillator |T]) 
at energy £ = 1 relevant for this work. The corresponding 
equipotential line has also been drawn in dashed line. 



in the calculations of this paper. They have been chosen 
as the semiclassically most relevant ones, in the sense 
that they are short, sym metric an d not too unstable (see 
discussion below, in Sect. Ill A 2). 



III. METHOD 

In this Section we describe the method that is used in 
our calculations. This description is made in three steps. 
First, we define the scar functions that are used to con- 
struct our basis sets. Second, we describe the procedure 
by which the elements of the basis set are selected and 
computed. And third, we introduce the mathematical 
tools that will be used to analyze the characteristics and 
quality of our results. 
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A. The scar functions 

This subsection consists of four parts. We first intro- 
duce (auxiliary) tube functions. Then some attention is 
paid to the details of the BS quantization procedure on 
the PO. In the third part, we discuss the actual calcula- 
tions of the scar functions, which are obtained by dynam- 
ically improving the tube functions, and constitute the 
primary ingredient of our basis elements. Finally, some 
examples of scar functions are presented and discussed. 

1. The tube functions 
Some auxiliary tube functions are first defined as 

C be (*,y)= Fdte^WtKxMt), (4) 
Jo 

where T is the period of the PO, and S n is the asso- 
ciated BS quantized energy (see next subsection). The 
wave function <fr(x, y, t) is a suitable wave packet, whose 
dynamics is forced to stay on the neighborhood of the 
scarring PO, (x t ,yt, P x ,t, Py,t)- For this purpose, we use 
a frozen Gaussian [33l|34] centered on the trajectory, that 
can be expressed as 

<j>(x,y,t) = exp {-a x (x - x t ) 2 - a y (y -y t ) 2 + 

(i/h)[P x , t (x - x t ) + P yit (y - y t )} + i lt ] , (5) 

where a x and a y define the widths along the two axis. 
In our calculation we take, for simplicity, a x = a y = 1, 
which is an adequate value for the problem that we are 
considering here. The phase 

lt = l -j\r(Pl T + Pl T ) - \» t 



is the difference between the dynamical phase of the or- 
bit [cf. Eq. (|3|] and a topological phase, proportional to 
the function which can be calculated by applying the 
Floquet's theorem [35]. For this purpose, the transver- 
sal motion is decomposed in the vicinity of a PO in two 
terms: one purely hyperbolic, describing the dilatation- 
contraction that takes place along the directions of the 
associated invariant manifolds, and another one, periodic 
in time, which is described by the matrix F(t) introduced 
in Eq. (11) of Ref. 23. After one period of time, li t = li 
is given by the winding number of the PO, which equals 
the number of half turns made by the manifold directions 
as they move along the PO. Also, this number is equal 
to the number of self-conjugated points plus the number 
of turning points existing on the PO. 

In general, ji t can be calculated as follows. First, the 
transversal monodromy matrix of the PO is computed 
and diagonalized. Recall that this matrix describes the 
linearized motion in the vicinity of a PO, and its eigen- 
vectors, £ s and £ n , give the directions of the stable and 



unstable invariant manifolds in that region. A trajectory 
starting on the stable/unstable manifold then approxi- 
mates/separates from the PO during the time evolution. 
The evolution of the eigenvectors £ s (0) and £ n (0) can be 
written, in the linear approximation, as follows 

Ut) = e~ xt F(t)U0), Ut) = e xt F(t)U0), 

where A is the stability index of the PO. After one pe- 
riod of time, £ S (T) and £ U (T) are e XT times shorter and 
larger, respectively, than the original vectors £ s (0) and 
£u(0), being the factor e XT equal to the absolute value 
of the largest monodromy matrix eigenvalue. Moreover, 
£ s (T) and £ M (T) are either parallel or antiparallel to £ s (0) 
and £ n (0), depending on whether the value of \i is ei- 
ther even or odd, respectively, i.e. the eigenvalues of 
the monodromy matrix are positive or negative. If they 
are positive (negative), the motion in the neighborhood 
of the PO is hyperbolic (hyperbolic with reflection) and 
F(T) is equal to the unit matrix, / (— I). Once the eigen- 
vectors, £ 3 (t) and have been calculated, we follow 
the evolution of the new vectors 

6(*) = e At &(*), Ut) = e~ xt Ut), 

which describe an unconventional motion in the vicinity 
of an unstable PO without hyperbolicity. For instance, 
let zpo(t) be the points of the PO as a function of time. 
Then, the evolution of a neighbor point z(0) = £po(0) + 
c 3 £ 3 (0) + c u £ u (0), with |c s |, \c u \ < 1, is given by z(t) = 
zpo(t) +c s £,s(t) + c n^n(^)- Finally, \i t can be obtained by 
following the angle swept by any of the previous vectors. 
Let us remark, that the value of \i t is not canonically 
invariant, in contrast with ii. 

2. The Bohr-Sommerfeld quantization rules 

The smoothing process implicit in the integration on 
Eq. Q renders a function with the probability density 
well localized along the PO. This localization effect over 
the PO is maximized, by constructive interference, when 
7t returns to the initial point with an accumulated phase 
that is a multiple of 2tt. This happens when the phase 
fulfills the so called BS quantization rule 

-^-^ = 27™, n = 0,1,2,... (7) 

At this point it is necessary to discuss the symmetry 
of the computed wave functions. The quantum eigen- 
states of the quartic oscillator ([!]) are classified according 
to the C^ v symmetry group, which has five irreducible 
representations (IRs), four of which (Ai, Bi, A^, B2) are 
one dimensional, and the other one (E) two dimensional. 
An elegant way to deal with this problem is to refer ev- 
erything to the fundamental domain defining the poten- 
tial. In our case this domain consists of the 1/8 region 
bounded by one semiaxis and the neighbor semidiago- 
nal in the case of the one dimensional representations, 
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A, £>, and the 1/4 region between the two semiaxis in the 
case of the two dimensional one, E. To translate this 
into semiclassical arguments, the POs must be desym- 
metrized by "folding" the original trajectories into the 
fundamental domain. In this way, eigenfunction symme- 
try characteristics turn into boundary conditions, Dirich- 
let (ip = 0) or Neumann (d±ijj = 0), at the axis (x, y = 0) 
and the diagonals (x = When dealing with POs, 

this is equivalent to introducing "artificial" hard walls 
boundaries in both the axis and the diagonals, which has 
two effects. First, they reduce the length, and then the 
topological (without the contributions arising from the 
desymmetrization) and mechanical actions in Eq. (|6| in 
an integer factor of p, given by the ratio between the pe- 
riods of the full and desymmetrized POs. Second, they 
have an additional more complicated effect in the Maslov 
index, which is different for the Dirichlet and Neumann 
cases, that has to be carefully taken into account. Ac- 
cordingly, the quantization condition ^ should be mod- 
ified in order to quantize a desymmetrized PO of period 
T/p, taking into account the appropriate boundary con- 
ditions of the PO. The new BS quantization rule then 
reads 



S T /P 

h 



[N D - N N ] - - \fi/p + N r }-= 27rn (8) 



where St/p, l-i/p + N r and N r are, respectively, the ac- 
tion, winding number, and number of reflections of the 
desymmetrized PO. Also, the number of excitations, n, 
has been 'reduced' to the fundamental domain. The num- 
ber of Dirichlet and Neumann conditions on the wave 
functions at symmetry lines (axis and diagonals) are 
given by and TVn, respectively. Obviously, N r = 
TVd + Nn. Furthermore, it can be shown that \i + 2pN& 
equals the Maslov index appearing in GTF [22j [36j [37] . 
A full discussion on the derivation of Eq. ( 8| can be found 
in Ref. [23j Finally, the semiclassically allowed BS quan- 
tized energies can be obtained by transforming Eq. ([sj) 
with the aid of the scaling relation ([3]), thus rendering 



£ 3 J A = 



2nh 

TsJb) 



4 2 



(9) 



where S = St is the action of the complete PO at en- 
ergy £ = 1. 

In Table [I] we summarize all the relevant dynamical 
information for the POs of Fig. [2] at £ = 1 [recall that 
they can be transformed to any other value of the energy 
by using the scaling relations in Eqs. ([2| and In 
the last column of the Table, we include an adimensional 
parameter defined as 



U = XTN a N u 



(10) 



measuring the relative relevance of each PO, in the sense 
that shorter, simpler, and less unstable orbits have lower 
values of 7Z. The integers N s and N t take into account 
the spatial and time-reversal symmetries of the orbits: 
N s corresponds to the number of different POs that are 
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Table I. Classical action S, stability index A, winding number 
/i, spatial and time reversal numbers N s ,N t , and relevance, 7Z, 
measured with Eq. ( |10[ ), for the POs of the quartic oscillator 
(l) shown in Fig. [2] at 8 = 1. 



obtained by application of the C± v symmetry operations, 
while Nt is equal to 1 when the PO is time-reversal and 
2 otherwise. Notice that the product N s N t equals the 
number of repetitions of the PO appearing in the sum- 
mation of the GTF, i.e. the number of similar POs that 
can be constructed at the same energy, which depends 
strongly on how symmetrical the PO is as well as on the 
IR that is being considered. 

Let us consider now the effect of desymmetrization. 
For the one dimensional IRs, the system is desym- 
metrized by reducing the configuration space to the re- 
gion x > y > 0. The corresponding desymmetrized POs 
are shown in Fig. [3j For the two dimensional representa- 
tion, the desymmetrized configuration space corresponds 
to the region x, y > 0, and the associated POs are plot- 
ted in figure Fig. |1J The corresponding information: p, 
TVd , and TVn , is given in Tables [TT] and |III[ respectively. 
In Table [TT] the reflections at the axis x = and the diag- 
onal x = y are considered, in addition to a reflection at 



the y axis for orbit 1. In Table III the reflections at the x 
and y axis are considered. We have separated the data in 
Ei and E^ components, corresponding to the cases sym- 
metric with respect to the axis x or y and antisymmetric 
with respect to the axis y or respectively. As POs 
numbers 2, 3, 10, 12, 14, 17 and 18 arrive at the origin 
x = y = forming a non-zero angle (see Figs. [3]and[4|, 
it is necessary in such situation to include reflections at 
the axes x and y simultaneously. 

Once the BS energies are calculated with the aid of 
the desymmetrized condition (pi), the corresponding tube 
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Figure 3. Same as Fig. [2] for the desymmetrized POs in the 
fundamental domain associated to the one dimensional irre- 
ducible representations. 

functions over the full PO (given in Fig. [5| are computed 
using Eq. @. These wave functions have n nodes in 
the desymmetrized region of the potential and Nb nodes 
along the symmetry lines, and are real only if the PO 
shows time-reversal symmetry. If this is not the case, a 
real function can always be obtained by combination of 
the tube functions obtained with two gaussian wave pack- 
ets given by Eq. ([5| running clock and counterclockwise 
along the orbit, respectively. Notice, however, that the 
time-reversal wave functions do not belong, in general, to 
any of the IRs of the system. Again, this represents no 
problem since proper symmetry-class adapted tube func- 
tions can be constructed combining the different wave 
functions obtained when the elements of the C± v sym- 
metry group act on the previously denned (real) tube 
functions. 



3. Computation of the scar functions 

The scar functions are finally computed as an improved 
version of the auxiliary tube functions introduced before 
that includes dynamical information of the system up to 
Ehrenfest time, T#. By restricting the integration time 
in this way, the geometrical support of the computed 
scar function is not only restricted to the PO itself but 
it also includes small pieces of the associated stable and 
unstable manifolds attached to the PO; see, for example, 
discussion in Ref. |24j 
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Table II. Period ratio p, and number of Dirichlet, Nb, and 
Neumann, TVn, boundary conditions relevant for the Bohr- 
Sommerfeld quantization (|9| of the desymmetrized POs in 
Fig. [3] in the one dimensional irreducible representations. 
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Table III. Same as Table |H| for the desymmetrized POs asso- 
ciated to the two dimensional irreducible representation, E, 
shown in Fig. [4] The data have been separated in two compo- 
nents, Ei and E2, corresponding to the cases symmetric with 
respect to x or y and antisymmetric with respect to y or x, 
respectively. 



The scar functions are calculated by propagating the 
corresponding tube functions Q, followed by a subse- 
quent Fourier transform for a finite lapse of time 

ili s n car (x,y) = J E dt cos (gp") g-ifA-fnlt/^tube^^^ 

(11) 

Te corresponds to the lapse of time that it takes to a 
typical Gaussian wave packet to spread over the area, A tr , 
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Figure 4. Same as Fig. [3] for the two dimensional irreducible 
representation, E. 



in a characteristic Poincare SOS of the desymmetrized 
system, and it is related to the Lyapunov exponent, A, in 
the following way 



2A V n 



(12) 



In our case A tr ~ 5.5278£ 3 / 4 and - 11.05555 3/4 for the 
one dimensional and two dimensional IRs, respectively, 
and A ~ 0.38485 1 / 4 , which does not dependent on the 
IR since it is associated to a generic chaotic trajectory of 
the system. 



The cosine window in expression (11) is intro- 
duced in order to minimize the energy dispersion, 

a = (^scar|^ _ £^2|^scar^ Q f the gcar f unct i ns [38], 

that can be semiclassically approximated as [27] 



7r h\(s 2 + AT E ) 

2 (5i + AT E )(52 + AT E )+5|' 



(13) 




Figure 5. (Color online) Probability density corresponding to 
the scar functions localized over the PO number 5 of Fig. [2] 
with n = 2, for the different irreducible representations of the 
system. The plot has been scaled to S = 1. The scarring PO 
is plotted superimposed. 



and then on A tr . This procedure effectively reduces the 
energy dispersion of the scar function with respect to the 
corresponding tube wave function by a factor that is pro- 
portional to \n(A tr /h). 



4- Some examples of scar functions 

Let us present now some representative examples of 
the scar functions that are obtained for the POs in Fig. [I] 
The value % = 1 will be used in all quantum computa- 
tions. Wavelets provide an efficient method to perform 
the time evolution appearing in Eq. ( 11 ), with a precision 
of at least six decimal places [46] . 



with si w 1.06078 and s 2 = tt/V2 - s t « 1.16066. No- 
tice that this magnitude depends on the IR through T#, 



In Fig. [5] we show the probability density correspond- 
ing to the scar functions constructed over the PO num- 
ber 5 with an excitation number n = 2, for all the IRs 
of the system. Since the guiding PO does not exhibit 
time-reversal symmetry, time-symmetrized frozen Gaus- 
sian functions, constructed combining wave packets ^ 
running clock and counterclockwise, have been used in 
the computations of Eq. Q. No additional spatial- 
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symmetrization is required in order to enforce the proper 
symmetry properties, since the time-symmetrized wave 
functions defined in this way are invariant under the ac- 
tion 

of the C^ v group operations. The BS quantized en- 
ergies also necessary in the calculations have been ob- 
tained with the aid of Eq. ([9| using the values given in 
Tables As it can be seen in the figure, the number 
of nodes of each of the wave functions is different, and 
equal to pN t (n + Nb/2). Furthermore, the functions be- 
longing to the one dimensional IRs are either symmetric 
or antisymmetric with respect to both the axes and the 
diagonals, those belonging to the two dimensional IR are 
on the other hand only symmetrical with respect to one 
axis and antisymmetrical with respect to the other. In 
panel (a), the A\ scar function having 16 nodes and lo- 
cal maxima over the axis and the diagonals (Neumann 
boundary conditions) is shown. On the other extreme, 
the A2 function [panel (b)] has nodes localized over these 
lines, exhibiting a total number of nodes of 24, 16 of 
which are "due" to the number of excitations, 4 come 
from the Dirichlet conditions at the axes, and the re- 
maining 4 nodes arises from the Dirichlet conditions at 
the diagonals. Similarly, the Bi (B2) scar functions have 
20 nodes over the PO, 4 of them due to the Dirichlet 
conditions at the diagonals (axes). Finally, for the E 
symmetry we have one Neumann condition over one axis 
and one Dirichlet condition over the other, this producing 
a total number of nodes of 10. 

In Fig. [6] we show the A\ symmetry scar functions 
constructed using the PO number 6 of Fig. [2] for dif- 
ferent values of the excitation number, n. As discussed 
before, these functions have been constructed imposing 
Neumann conditions on the boundaries of the desym- 
metrized POs, so that n is equal to the number of nodes 
in the fundamental domain, i.e. 8n in total. The local- 
ization on the scarred PO as n increases is notorious. 

In Fig. [7] we present some results for the energy disper- 
sion of our scar functions, together with the correspond- 
ing semiclassical estimation obtained from Eq. ( fl3| ) . As 
can be seen, our functions are very well localized in en- 
ergy, and their dispersion grows moderately with it. This 
is a key point for the aim of this paper, since it allows 
the definition of a very efficient basis set. 



B. Definition of the scar functions basis set: 
Selective Gram-Schmidt method 

The second pillar of our method is the definition of 
the selection procedure of the scar functions forming the 
basis set, that is subsequently used in a standard diago- 
nalization of the associated Hamiltonian matrix to obtain 
the eigenstates of the system. 

To define our basis set, we have generalized the usual 
Gram-Schmidt method (GSM) [39 , and developed a new 
selective Gram-Schmidt method (SGSM) able to choose 
a basis set of linearly independent functions in a vecto- 
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Figure 6. (Color online) Scar wave functions with A\ symme- 
try constructed using the PO number 6 of Fig. [2] for different 
values of the excitation number, n. The scarring PO has been 
superimposed. 



rial space from a larger (overcomplete) set of functions, 
that can be used to efficiently compute chaotic eigen- 
functions of our system in a given energy window. The 
SGSM is specially useful when computing highly excited 
eigenfunctions in a small energy window, since the size 
of the basis set reduces considerably in this case. The 
procedure starts from an initial set of N scar functions, 
|^ (0) ), from which the SGSM selects the minimum num- 
ber of them, Nb < TV, necessary to adequately describe 
the Hilbert space defined by the eigenfunctions whose en- 
ergies are contained in the energy window, i.e. the SGSM 
defines a basis set in that window. The elements of the 
basis set where subindex i orders the elements 

according to their semiclassical relevance (see discussion 
below), are automatically selected with the aid of the 
conventional GSM. Thus, associated to the basis 
we will construct an auxiliary basis \(fi), formed by the 
orthogonalization of For example, if we set 

\<pi) = i^f) 
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Figure 7. (Color online) Dispersion of the scar functions of symmetry A\ (red plus signs), A2 (green crosses), B\ (dark blue 
stars), B2 (pink open squares), and E (light blue filled squares) constructed over the periodic orbits 1-6 shown in Fig. [2] The 
results obtained with the semiclassical approximation ( 13 ) are plotted superimposed in black continuous line. 



then, the second auxiliary function \(f 2 ) is given by 



1^ = T^T 



(l)l 



where 22 7^ ji and 
and so on. 

In our SGSM method, the selection procedure of the 
basis functions of a given symmetry for the calculation 
of the eigenenergies, £ , contained in the energy window 
defined by 



£~ <£< £+ 



(14) 



is done automatically by using a definite set of rules, 
which are based on a semiclassical selection parameter, rj. 
This parameter is defined in such a way that it takes 
into account in a simple form the dispersion of the scar 
functions, the simplicity of the PO, and the density of 
states of the system (which is only relevant when the 
energy window is large). For a given scar function, rjj is 
given by 



(SSjf^TjNsjN^ 



j ' 



(15) 



where pj is the mean density of (symmetry-class) states 
at the BS energy of interest (£j), aj is the dispersion of 
the scar function, and S£j is defined as 



5£i 



£~ 
0, 



£ ) , if £j < £ 
iff 



< £j < £+ 



(16) 



£j -£+ if£j >£+. 



Other criteria could be used to define the parameters in- 
troduced in Eqs. (10) and (15). For example, one could 



drop out the function S£ ) appearing in Eq. ( 15 ) and still 



get quite accurate results, specially in large energy win- 
dows. However, this function is included to improve the 
numerical accuracy by reducing the boundary effects. We 
thus believe that the previous equations are very straight- 
forward in order to substitute the contribution of the 
longer POs appearing in the GTF by the interaction of 
the shortest ones, following the short PO theory devel- 
oped by Vergini et al. [22j |23l [27] , although other defini- 
tions are possible. 

The SGSM is then defined in an algorithmic way as 
follows: 



O.a With the method described in Sect. |III A[ we 
compute all normalized scar functions, |^°^) with 
the smallest values of 1Z and BS quantized energies 
£j in the enlarged energy window 



£- 



2a < £, < £^ 



2a, 



(17) 



where a is given by Eq. (13) for A = A. This is the 



most time demanding step of the procedure. 

It can be a priori expected that the overlap of the 
scar functions outside this enlarged window with 
the desired system eigenfunctions is negligible, due 
to the fact that they were constructed minimizing 
their energy dispersion. 

• O.b The semiclassical selection parameter, rjj, of 
Eq. (15), associated to each scar function in this 



first approach to the final basis set, is then calcu- 
lated. 

• 1. From this initial set of scar functions, 1^°^), we 
select a smaller number of them, < N, forming 
the basis set that is optimal for our purposes, in 
the following way. Notice that the number of scar 
functions calculated in this way, TV, should always 
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be greater or equal to 

N b = N sc (£ + + 2a) - N SC (S~ - 2a) + c b ap, (18) 

where, N sc (£) is the semiclassical approximation 
to the number of states with an energy smaller 
than £ , and the term c b ap enlarges the window 
size to take into account border effects. If this is 
not the case, more (longer) POs, and consequently 
more scar functions, must be included at this step. 

The first element of our basis set is the scar function 
with the smallest rjj value 



l*S?>> 



with 



1 



(19) 



This choice gives priority to the scar functions 
which are more localized in energy over simpler 
POs, i.e. with shorter period and being more sym- 
metric (smaller N S j and N t j). In a similar way, we 
will denote from now on by \cpj) the auxiliary func- 
tions necessary for the selection of the scar func- 
tions. 

2. a The remaining scar functions are then orthog- 
onalized to function IV^) using the usual GSM 



• n.b The n-th basis element, is then selected 

as the one for which 



O-!) |2 



IV; 



(™-l)|2 ' 



Vj 



(24) 



j 9^ J 1:32,- .., 3 ] n 



and then the next auxiliary function is constructed 
as 



\<p n ) = 



i^r 1} > 



|V; 



(n-l). 



(25) 



• The procedure finishes when the number of se- 
lected elements of the basis set equals N b given by 
Eq. pi. 



Finally, the corresponding Hamiltonian matrix is com- 
puted in the basis set of scar functions, or alternatively 
in the equivalent basis set of auxiliary functions, with the 
help of the wavelets, which provide an accuracy of at least 
14 decimal places for the matrix elements. The diagonal- 
ization using standard routines [40] gives N b eigenstates 
in the energy window defined in (14). 



ivf> 



h&i 0) > 



<^ 7 (0) >K>, 



3*31- (20) 



2.b The second element of the basis set is \ip 
where the index 22 satisfies 22 7^ ji •> and then 



(0)\ 



Vj2 



= max 



Vj 



(21) 



37^31 



where the norm in the numerator has been intro- 
duced in order to make the basis set elements as 
different as possible between them. Indeed, notice 
that after the orthogonalization condition (20) the 



more similar is to \(fi), the smaller the 



norm of the function • Then the auxiliary 

function ^2) is computed as 



^2) 



1^1 



(22) 



The previous steps, 2. a and 2.b, are repeated for all 
the remaining basis elements in the initial basis set 
of scar functions, in such a way that the n th step 
in the procedure is defined as: 

• n.a New functions are obtained by orthogonal- 
ization to the auxiliary one in the previous 
step, |^ n _i), 



1^ 



(n-l) x 



- 2 )\ 



3 + jl 



(n-2) x 



32, 



l^n-l), 
••i 3n—l- 



(23) 



C. 



Local representation, scar intensities, and 
participation ratio 



In order to analyze the performance of our method we 
will use in this work a local representation, |^ oc ), defined 
as 



\N) 



C Nj\¥ l j OC ), 



(26) 



where the coefficients C^j = ((p l j OC \N), to reconstruct the 
eigenfunctions of the system in such a way that we can re- 
gain the attractive initial intuitive interpretation. Again, 
the procedure to compute the different |^ oc ) functions 
will be described here in an algorithmic way. 



• 1. The first element of the local representation 
is taken as the scar function, that has the 

largest value of the scar intensity, x±, defined as 



» 



|(^ n) |7V)| 2 . 



(27) 



That is 

loc\ _ 1 /,(0) 



l^° c ) 



|^}/), with x- 



(0) 



X 3l 



max{x 0) }, (28) 



• 2. a In order to identify the second largest scar in- 
tensity, one has to calculate the orthogonal part 
of the remaining scar functions 1^°^) to |^i° c ) by 



computing 



ivf> 



h&i 0) > 



(^'riVf )ki° c ), 3*3i- (29) 



Joe 
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• 2.b The second element of the local representation 
is taken as 



Ml 

31 1 



(30) 



with X2 



JX) 



The procedure is then continued, so that the n-th 
element of the local representation is computed in 
a similar way: 



n.a The orthogonal part 

l^ n_2) >iWiA,...Jn-i tO \<Pn-l 

computing 



of the functions 
) is calculated by 



\4 n ~ 1} ) = l^"" 2) ) - (^°-il^ n - 2) )l^°-i), (3D 



• n.b and the n-th element of the local representation 
is given by 



i^-r 1 ') 



(n-l). 



(32) 



with x n = x 



(n-l) 



max{x f \j ^ ji, j 2 ,-.., jn}. 



Let us remark here that the scar intensities, Xj, pro- 
vide information on the localization properties of the 
eigenfunctions. However, although the largest scar inten- 
sity, xi, provides a faithful information on the localiza- 
tion on l^j^), the same is not true of \^f^) since in their 
construction, the contribution in the subspace spanned 
the functions defined previously in the SGSM was sub- 
tracted. Nevertheless, the sum x\ +X2 + ... + x n provides 
information about the projection of \N) onto the sub- 
space defined by |^), l^), • • • , l^)- 

Let us finally present some useful results concerning 
the scar intensities defined in Eq. (27). Assuming that 



the distribution of these magnitudes follows a Gaussian 
law for chaotic eigenfunctions, a semiclassical approxi- 
mation for the average can be obtained, as discussed in 
Ref. Hi] Indeed, the averaged value of the j-th largest 
scar intensity is given by 



H 1 

0~ r \ 7T 



1/2 



*3 ~ ln a 3 
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2 



, (33) 



where 



o~ r per 7v , ctj = Zj 



In 



y/2a r _ \n(aj + 287/128) 



J 



17/8 



(34) 

being z a random variable with averages Z\ c± 0.577 
and Z2 — 13/48 for the two largest scar intensities. Let 
us remark that Xj goes, respectively, to zero and infin- 
ity for large and low values of the energy. Also, the 
above defined semiclassical expression can be improved 
by including a higher order correction term in <x,, so 



that aj = Zj + \n(y/2a r /j + Cj). Adequate values for Cj 
are obtained by fitting to actual quantum calculations 
of Xj. Notice that these corrections are of order 0(l/oy), 
while ctj = 0(1/ \na r ). 

We conclude this subsection by considering the par- 
ticipation ratio, i?, that is defined, taking into account 
expression (26), as 



R 



3 = 1 ^Nj 
3=1 ^Nj 



(35) 



This magnitude gives an idea of the number of basis el- 
ements that are approximately necessary to reconstruct 
the original eigenstate, \N). Accordingly, this is a pa- 
rameter very relevant for our discussions, since it can be 
used to compare the quality of two different basis sets. 
Namely, the lower the value of the R, the better the basis. 



IV. RESULTS 

In this section we present and analyze our results on 
the use of scar functions as basis sets for the calculation 
of eigenstates of classically chaotic systems. 



A. Calculation of the eigenstates 

The lowest eigenvalues and eigenfunctions of the 
Hamiltonian operator corresponding to the quartic os- 
cillator defined in Eq. have been calculated using a 
basis set of scar functions constructed using the SGSM 
defined above. 

For this purpose, the 18 POs presented in Fig. [2j cho- 
sen using the relevance parameter defined by Eq. ( |To| ), 
have been used. From them, an initial basis set of scar 
functions was constructed, defined by the reference en- 
ergies £~ = and £+=135, 150, 140, 140, and 82, 
with a=0.5433, 0.5526, 0.5465, 0.5465, and 0.4622, for 
the symmetry classes Ai, A<i, Bi, and E, respec- 
tively. This initial basis set consists of ~900 scar func- 
tions for each of the one dimensional symmetry classes: 
Ai, A 2 , £?i, £?2, and also for E\ and E 2 . The associated 
energies were obtained from the BS quantization rule ([9| 
using the parameters given in Tables [TpTT| The lowest 
values corresponding to the A\ symmetry class, rescaled 
with a power of 3/4 so that they appear equally spaced, 
are represented with thin black lines in Fig. [8] The next 
step in our procedure is the reduction of this initial basis 
set. This is done by applying the SGSM described in the 
Sect. |IIIB| using a value of c b = 2 in Eq. ([18]) . With this 



criterion, ~420 scar functions of each one dimensional 
and Ei and E% symmetry classes are automatically se- 
lected, out of the whole set. The energies corresponding 
to these selected A\ scar functions have been highlighted 
with thick blue lines in Fig. [§J 

Using the resulting final basis set, we have computed 
the corresponding Hamiltonian matrix and, by direct di- 
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agonalization, the eigenenergies and eigenfunctions of the 
system. The results for the 270 low-lying states of A\ 
symmetry are shown in the leftmost (red lines) tier of 
Fig. [8j In Table IV the 50 low-lying numerical values 
for all symmetries are reported. Several comments are 
in order. First, a careful comparison reveals that the 
first ~400 of each symmetry have the same accuracy (not 
only the eigenenergies but also the eigenfunctions) that 
is obtained with other standard methods pQ, and then 
the number of well converged states in our calculation 
is of the same order as the number of elements in the 
basis set. Second, some first qualitative conclusions re- 
garding the excellent performance of our basis set can be 
obtained from a careful consideration of the results pre- 
sented in Fig. [8] To this end, recall again the extremely 
low dispersion of our scar functions, which means that 
they minimally spread among (or contribute to) states in 
the eigenenergy spectrum far from their BS quantized en- 
ergies. For example, for the most excited states only ~6 
scar basis functions are needed for a satisfactory descrip- 
tion. This number dramatically decreases for smaller ex- 
citations, getting as low as ~1 for states in the interval 
AS ~ — 5, ~2 for states in the interval AS ~ 5 — 10, 
and so on as the scar states get "bright", being high- 
lighted in blue in our plot. This argument will be made 
more quantitative in the rest of the discussions presented 
below in this section, and particularly when the eigen- 
state participation ratios are considered. Finally, it is 
worth emphasizing that besides the type of calculation 
presented here, namely, the computation of the low- lying 
eigenstates, the SGSM introduced in this work can be ad- 
vantageously used in the computation of eigenstates in an 
energy window, i.e. £~ ^ 0, something which is specially 
useful to compute only highly excited eigenfunctions with 
small basis sets. This is something that we have tested 
for different energy ranges. In the next subsection we 
present the results associated to the eigenfunctions | N)a 1 
with N = 271 — 282, calculated restricting the calcula- 
tion to the energy window 100 < £ < 103 with a basis set 
formed by only 25 scar functions. It is quite impressive 
that using such a small basis set the error in the energy is 
smaller than 0.12 in units of the mean level spacing (ac- 
tually, it is even smaller than 0.06 for all the computed 
eigenvalues except for |279)aJ, whereas the overlap of all 
the computed eigenfunctions with the corresponding ex- 
act ones exceeds 99.8%. Further details on these results 
will be presented in the following section. 



Figure 8. (Color online) Scaled energies for the A\ states 
of the quartic oscillator |l]): (red leftmost tier) numerical 
eigenenergies, (thin black lines) BS quantized energies, (thick 
blue lines) energies of the scar functions selected by the 
SGSM. 



B. Reconstruction of the eigenfunctions 

Let us now analyze the results obtained in the previous 
subsection. In the first place, we discuss the structure of 
some representative examples of the eigenfunctions ob- 
tained for the quartic oscillator ([!]) with our basis set of 
scar functions, by examining their reconstruction using 
the local representation described in subsection |III C| 

We start by the simplest case, which corresponds to 
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N A! A 2 Bi B 2 E 

1 0.56323 4.1023 1.6175 2.5230 1.2241 

2 1.8848 5.9503 2.7455 4.7256 2.2570 

3 2.8638 7.5442 3.8675 6.0561 3.2537 

4 3.8563 8.9638 5.0139 7.2567 3.6376 

5 4.8286 9.7951 6.1773 8.2227 4.4506 

6 5.2584 10.658 6.8364 9.1029 5.1299 

7 6.2126 11.985 7.4816 10.468 5.5775 

8 7.4115 12.861 8.7094 11.352 6.2642 

9 7.9052 13.616 9.3168 11.914 6.8080 

10 8.6947 14.793 10.080 12.708 6.9510 

11 9.3055 15.393 11.188 13.637 7.9451 

12 10.087 16.069 11.534 14.276 8.1741 

13 10.664 16.731 12.525 15.134 8.3392 

14 11.452 17.702 12.908 15.928 9.3842 

15 11.960 18.287 13.550 16.571 9.4158 

16 12.790 19.115 14.309 17.215 9.8674 

17 13.140 20.011 14.934 18.242 10.448 

18 14.298 20.467 15.827 18.702 10.804 

19 14.714 21.183 16.234 18.976 11.013 

20 15.003 21.789 17.007 19.988 11.520 

21 15.831 22.345 17.391 20.220 12.004 

22 16.024 23.217 18.099 21.270 12.267 

23 16.924 23.673 18.928 21.821 12.786 

24 17.384 24.062 19.068 22.073 13.091 

25 17.989 24.810 19.885 23.038 13.526 

26 18.517 25.348 20.202 23.556 13.814 

27 18.972 26.495 20.604 24.428 14.109 

28 19.905 26.629 21.668 24.591 14.541 

29 20.184 27.083 22.118 25.140 14.729 

30 20.592 27.860 22.358 25.967 15.069 

31 21.163 28.462 23.087 26.104 15.646 

32 21.931 28.946 23.754 26.682 15.893 

33 22.228 29.614 24.066 27.357 16.209 

34 22.458 29.711 24.939 27.950 16.600 

35 23.247 30.334 25.079 28.478 16.801 

36 23.727 31.099 25.513 29.093 17.373 

37 24.010 31.516 25.916 29.241 17.407 

38 24.653 32.199 26.817 30.233 17.768 

39 25.231 32.626 27.121 30.333 18.152 

40 25.576 32.831 27.443 31.015 18.339 

41 26.121 33.681 28.180 31.511 18.764 

42 26.530 34.089 28.491 32.184 19.147 

43 27.098 34.516 28.932 32.312 19.174 

44 27.403 35.335 29.345 33.214 19.765 

45 28.097 35.928 29.980 33.441 20.005 

46 28.727 36.002 30.711 34.020 20.136 

47 28.946 36.551 31.076 34.323 20.642 

48 29.342 36.933 31.467 34.651 20.898 

49 29.475 37.373 31.568 34.953 21.120 

50 30.225 37.855 32.483 35.880 21.440 
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Figure 9. (Color online) (Right tier) Probability density for 
the eigenfunctions |208) Al , |291) Al , and |303) Al obtained 
from our variational calculation using a basis set of scar func- 
tions, and (left tier) the same functions reconstructed using 
only one scar function |PO,n) Al of the basis set. The scar- 
ring PO has been plotted superimposed in panels (a)-(c). The 
energy is shown on the left corner of each panel, and the over- 
laps of the eigenfunctions with the scar functions (left tier) 
and the participation ratios (right tier) are given on the right 
corner of the panels. 



states with eigenfunctions that appear strongly scarred in 
the sense discussed by Heller in Ref.|6j This happens, for 
example, for states |208) Al , |291) Al , and |303) Al , which 
are highly localized over POs numbers 3, 8, and 1 of 
Fig. [2j respectively. The probability densities are shown 
in the right tier of Fig. [9j In all cases, the first, and al- 
most exclusively contributing, element of the reconstruc- 
tion is a scar function corresponding to the PO scarring 
the eigenfunction, as could be expected a priori. These 
scar functions, which will be labelled |PO,n) x with PO 
indicating the number of the orbit in Fig.|2j n the number 
of excitations along it corresponding to the BS quantiza- 
tion condition (|9|, and \ the IR, are presented in the 
left tier of the figure. The PO has been plotted super- 
imposed to the corresponding probability density. The 
associated energies are given in the lower left corner of 
each panel. In the lower right corner of the panels (a)- 



Table IV. Eigenenergies, £, for the eigenstates of the quartic 
oscillator (HI obtained with our basis set of scar functions. 
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Figure 10. (Color online) Same as in Fig. [9] for the eigenfunc- 
tion |243)ai i 3 )^ which requires two scar functions, |9, 50)ai 
(b) and |6, 29) a x (c) for its reconstruction, which is shown in 
(d) and (e). The information included in the different panels 
is the same as in Fig. [9] 



(c), we have indicated the value of the overlap between 
the eigenstate and the basis set scar function. As can be 
seen, this overlap is always larger than 86.0% in all cases 
considered here. Furthermore, the fact that these states 
are well represented by only one state of the scar basis 
set makes the values of the corresponding participation 



ratio very small, as discussed in the subsection |IVD| The 
values of these participation ratios are given in the right 
corners of panels (d)-(f) in Fig. [9] 

Let us consider next eigenstates with more complex 
structure. This is the case, for example, of |243)ai shown 
at the top of Fig. [To] As can be seen in panel (a), 
the eigenfunction for this state is concentrated on a sin- 
gle PO, i.e. PO number 9, but the localization is not 
as strong as in the previous examples. Actually, the 
reconstruction of this eigenstate requires the combina- 
tion of at least two scar functions, as the results pre- 
sented in the other panels indicate. Indeed, the scar ba- 
sis function |9,50)a 15 shown in panel (b), only accounts 
for 75.0% of the eigenstate, while when function |6, 29) a x 
[cf. panel (c) of Fig. 10 is included, this figure raises 
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to an acceptable 92.6%~As discussed before, the corre- 



Figure 11. Same as Fig. 10 for |201)ai which requires four 
scar functions for its reconstruction. 



sponding value of the participation ratio is expected to 
be larger, R = 1.68, in this case. 

A similar example, but with an eigenfunction exhibit- 
ing an even more complicated structure is shown in 
Fig. 11 In this case, the eigenfunction |201)ai, which 
is localized over the "box" PO number 5 [see panel (a)], 
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requires the combination of at least four scar functions 
localized over different POs [see panels (b)-(e)] for an 
adequate reconstruction. The relevant figures for this re- 
construction, i.e. overlaps and participation ratios, are 
indicated in the figure, and the same comments made 
before apply to this case. 

To conclude this series of examples, let us consider 
finally the most general case consisting of (non-scarred) 
eigenstates showing the very irregular patterns character- 
istic of chaotic states 142 . For this purpose, we present 
in Figs. fl2] (a) and Il3] (a) the probability density for 
states \2T1)a 1 and |327)ai, respectively. The complex 
nodal structure inherent to these wavefunctions can be 
unfolded, however, by our SGSM which reveals the im- 
portance of each PO, thus proving a dynamically oriented 
analysis of them. As can be seen, the two cases consid- 
ered here are essentially reconstructed by combining only 
five scar functions, which are also shown in order of im- 
portance in the two figures [panels (b)-(f)]. 

From the data contained in the figures, notice how in 
all the examples presented here the scar functions giv- 
ing the largest contributions to the reconstruction of a 
given eigenstate are those whose BS quantized energies 
are closer that of the corresponding eigenenergy. 

We close this subsection by presenting in Table [V] the 
structure of the eigenfunctions calculated in the energy 
window 100 < £ < 103. These eigenfunctions are shown 
in Fig. [14] . As in the previous examples, most of the 
probability density is reconstructed by combining few 
basis elements, thus demonstrating the efficiency of our 
method for the calculation of excited states. 

Further results on the structure of the eigenfunctions 
of the quartic oscillator in our scar function basis set 
are systematically presented in the supplemental mate- 
rial 031. * 



C. Scar intensities 

A more global idea of how the eigenfunctions are re- 
constructed by the scar functions can be obtained by 
considering the corresponding scar intensities defined in 
Eq. (27). In the two panels of Fig. 15 we show (pink plus 
signsjthe variation with the energy of the two largest 
scar intensities, x\ and #2, for the computed A\ eigen- 
functions. Recall that these quantities give the (squared) 
contribution of the two most important scar functions in 
the local representation of each eigenfunction. As can 
be seen, both x\ and x^ fluctuate violently, thus appear- 
ing rather dispersed in the figure. To clearly observe 
the tendency in the data we have also plotted the corre- 
sponding mobile mean (blue crosses), computed using 20 
points. The results indicate that the largest intensity, xi, 
starts from a very high value, and then monotonically de- 
creases, first very rapidly up to £ ~ 22, and then much 
slower, never getting in mean below ~0.4 in the energy 
range that we are considering. Notice that the points 
where x\ is much larger that its average value corre- 




Figure 12. Same as Fig. 10 for |277) Al 
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Table V. Reconstruction of the eigenstates \N)a 1 cal cula ted in the energy window 100 < £ < 103 with the 25 scar functions 
|PO, ti)a 1 - The participation ratio, R, defined in Eq. (35) and the accumulated reconstruction overlap, £%, are also given. 



N R PO,n,S% 



271 


4.81 


9,53,40.5 


9,53,40.5 


4,64,56.8 


8,31,62.6 


18,45,69.5 


14,49,74.0 15,47,77.9 5,13,81.4 10,52,84.5 17,39,87.4 


272 


3.53 


15,47,48.4 


13,35,62.2 


9,53,76.7 


10,52,82.4 


4,64,86.5 




273 


2.09 


8,31,67.6 


6,31,76.6 


5,13,87.1 








274 


2.00 


1,54,69.7 


5,13,78.2 


8,31,84.0 


18,45,88.1 






275 


5.24 


2,54,31.0 


5,13,55.9 


8,31,69.0 


13,35,76.0 


4,64,84.9 


17,39,87.5 


276 


3.36 


6,31,28.7 


8,31,72.6 


10,53,85.5 
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4.68 


10,53,29.5 


9,54,51.5 


4,65,78.5 


5,13,83.8 


3,21,86.5 
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5.27 


10,53,30.8 


3,21,53.3 


4,65,71.4 


17,40,78.2 


18,46,82.0 


6,31,84.6 8,31,86.3 
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7.85 


9,54,21.8 


4,65,34.2 


15,48,40.1 


3,21,45.5 


13,36,52.5 


5,13,56.1 2,55,58.5 1,55,79.4 1,54,81.1 17,40,81.9 






18,46,88.1 
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2.08 


9,54,67.7 


15,48,71.3 


16,43,73.9 


10,53,76.5 


5,13,77.5 


18,46,78.3 13,36,79.2 2,55,79.9 1,55,92.6 


281 


3.30 


17,40,49.3 


15,48,69.7 


9,54,77.9 


4,65,85.4 






282 


4.20 


3,21,44.3 


13,36,52.7 


15,48,65.7 


11,52,73.7 


4,65,82.0 


10,53,84.1 9,54,86.1 



spond to scarred states in the sense of Heller [6 j. Si- 
multaneously, the second largest contribution, #2, first 
increases up to ~0.2 at £ ~ 22, and then decreases ex- 
tremely slowly (notice that although this behavior in not 
noticeable in the plot this is true since x\ > x^). More- 
over, the values given by the semiclassical approxima- 
tions (33), without and with the higher order correc- 
tion terms with c\ ~ — C2 — 0.30, have also been plot- 
ted superimposed in red solid and green dashed lines, 
respectively. The agreement between the quantum and 
semiclassical calculations, particularly in the second case 
(green dashed line) is rather good, except for low values 
of the energy due to the (unrealistic) singularity inherent 
to Eq. ([33]). 



Similar results for the A2 eigenstates are shown in 
Fig. 16 Here, only the values of the averages and the 
semiclassical estimates are shown for simplicity. The 
same comments made before for the A\ symmetry ap- 
ply- 



D. Participation ratios 

Also interesting for the analysis of the eigenstate struc- 
ture is the consideration of the participation ratios dis- 



cussed in subsection III C To give an idea of the expected 
bounds for these variables let us consider two extreme 
cases. First, the minimum value of the participation ra- 
tio is equal to R m i n — 1, which is obviously obtained 
when the system eigenfunctions are used. On the other 
extreme, a large value of the participation ratio corre- 
sponds to a situation when all basis functions signifi- 
cantly contribute to the description of each eigenstate. 
One such case, that retains however the character and 
then the efficiency of a semiclassical description, is that 
of Ref. 4. There, all calculations are done in a charac- 
teristic Poincare SOS of an ergodic system. Accordingly: 
1) the size of the basis set, A/5, can be estimated by the 



Weyl expression iV& = A tr /(27rh), being A tr the area in 
the SOS, and 2) all the coefficients in a normalized expan- 
sion are expected to behave as random complex indepen- 
dent variables following a Gaussian distribution of zero 

1/2 

mean and dispersion l/N b f . Then, it is straightforward 
to show that 



R« 



A tr 



(36) 



Let us remark that R sc is still a reasonably small value, 
since it is associated to an optimized, in the semiclassi- 
cal sense, basis set, thus giving participation ratio values 
much smaller than, for example, for other frequently used 
basis sets, such as harmonic oscillator products [1 , or a 
discrete variable representation [44] . 



In Fig. [TT| we plot with pink plus signs these partici- 
pation ratios in the local representation as a function of 
the energy for the A\ (top) and the A2 (bottom) symme- 
try eigenfunctions, respectively. Similarly to what hap- 
pens with the scar intensities we find a wildly oscillatory 
behavior. Accordingly, we consider the averaged value, 
computed again as a mobile mean, which is also shown in 
the graph using blue crosses. The results obtained with 
the semiclassical expression (36) have also been plotted 



for comparison with dashed green line. As can be seen, 
for energies up to £ ~ 22 the values of the participation 
ratios are very small, and always lower than 3. As en- 
ergy increases the mean participation ratio also increases, 
but very moderately, and it can be accurately fitted with 
the simple expression (also shown in the figure with red 
continuous line) 



R = (cr r 



(37) 
Here 



in terms of the variable a r defined in Eq. (34) 
C = 8/3 for the one dimensional symmetry representa- 
tions, A and and ( = 2 for the two dimensional one, 
E. The moderate increment behavior found is an indi- 
cation of the quality and good performance of the scar 
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Figure 13. Same as Fig. 



10 



for |327) Al . 



Figure 14. (Color online) Quartic oscillator eigenfunctions 
|iV)Ai corresponding to the reconstruction data of Table |v| 
Same coordinates and scaling as in Fig. [5] 
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Figure 15. (Color online) Largest scar intensities x\ (top) 
and X2 (bottom) (pink plus signs) for the eigenfunction with 
A± symmetry obtained in our calculation with a basis set 
of scar functions. The average value, computed as a mobile 
mean, is plotted with the superimposed with blue crosses. 
The semiclassical estimates obtained with Eq. (33) for c\ — 
C2 — and c\ c± — C2 — 0.30 are also plotted with red solid 
and green dashed lines, respectively. 



basis set we have used in our calculation, and shows how 
a given eigenstate can be reconstructed with just a few 
number of basis functions. Actually, it has been numer- 
ically checked that the cases with the highest values of 
the participation ratio, this number can be substantially 
reduced by increasing the size of the basis by including 
more POs. 

The participation ratio results can be further analyzed 
by considering the associated statistics. For this purpose 
we define the adimensional coefficient 



R- 1 



(38) 



which is a positive random variable that is found to follow 
the Weibull distribution [45l. 



AvW = f (f ) 



k-1 



-{x/l) k 



(39) 



k and / being fitting parameters. The associated accu- 
mulated distribution can be easily calculated as 



W w (r 



Jo 



Pw(%) dx = 1 — e 



-(r/l) k 



(40) 




Figure 16. (Color online) Average value of the largest scar in- 
tensities x\ (blue upper plus signs) and X2 (pink lower crosses) 
of the eigenfunctions of A2 symmetry. The semiclassical esti- 
mates obtained with Eq. (33) are also plotted with red solid 



and green dashed lines, respectively. 




Figure 17. (Color online) Participation ratio in the basis set of 
scar functions (pink crosses) for the eigenfunctions of A\ (top) 
and A2 symmetry, as a function of the energy. The average, 
computed as a mobile mean, is shown with blue crosses, and 
it is well approximated wit h the continuous red curve given by 
the simple expression (37). The dashed green line represents 
semiclassical estimate given by Eq. (36). 



The corresponding results for both the computed and 
fitted accumulated distributions W(r) and WwW for 
the eigenfunctions of A\ and A2 symmetry are shown 
in Fig. 18 The difference between W(r) and WwW, or 
error, is also shown in the inset. As can be seen, this 
difference is very small for both symmetries, as it is also 
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Figure 18. (Color online) Accumulated probability distribu- 
tion, W(r), of the participation ratio adimensional coefficient, 
r, for the eigenfunctions of A\ (red upper points) and A 2 
(green lower points) symmetry obtained in our calculation 
with a basis set of scar functions. The corresponding val- 
ues of the Weibull distribution, WW(r), for l Al = 2.8323 ± 
0.0003, k Al = 1.2030 ± 0.0002, l M = 3.1867 ± 0.0001, k A2 = 
1.1753 =b 0.0001 are shown in continuous blue and dashed pink 
lines, respectively. The difference W(r) — Ww(f) is shown in 
the inset. 



the case for the rest of the symmetry classes. 



E. Upper bounds to errors in eigenenergies and 
eigenfunctions 

In this subsection we obtain expressions for the upper 
bound to the error in the eigenenergy and eigenfunctions 
obtained in our calculation versus their dispersion [27] . 
These upper bounds are specially useful in the calcula- 
tion of highly excited states, when large basis sets are 
required. 

Usually, the convergence of approximate eigenstates, 
both in energy and wave functions, is assessed by com- 
paring the results obtained using two basis sets of differ- 
ent size. In our case, we will take as the 'exact' results, 
denoted by £' N , \N f ), those obtained variationally by di- 
agonalization in a very large basis set (~5000 elements) 
of harmonic oscillator eigenfunctions [TJ [2] . 

In Fig. [19] we present such errors for both eigenenergies 
(bottom panel) and eigenfunctions (top panel) of all sym- 
metry classes, computed as A£ r = \£n~ &n' \p (measured 
in mean energy level spacing units of each symmetry- 
class) and 1 — (N\N f ) 2 , respectively, as a function of the 
reduced dispersion of the eigenfunctions 

cr r a N p, 

where a at is the dispersion of the eigenfunction \N) X . As 
can be seen, both A£ r and 1 — (TV | TV') 2 follow to a great 
degree of accuracy a power law for <j r < 0.3, thus giving 
upper bounds that can be expressed as 




1 G~ 1 ~ i 1 1 1 1 ill i i 1 1 1 1 mI i i 1 1 1 1 mI 



j i i 1 1 1 in i i i 1 1 1 n 



0.001 0.01 
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Figure 19. (Color online) Error in the eigenenergies (bottom 
panel) and eigenfunctions (top panel) of the eigenstates ob- 
tained in our calculation using a basis set of scar functions, 
estimated as described in subsection |IVE| as a function of the 
relative dispersion, a r — (JnP- The continuous lines indicate 
the upper error bounds given in Eq. (41). 



There is a good reason for calculating the previous upper 
bounds as a function of <j r . Indeed, this latter magnitude 
can be computed straightforwardly, so that one can es- 
timate the errors in the results with the aid of Eq. (41), 



^2.5 



AS r < 



1 - (N\N'} 2 < 



^2.5 



100 



(41) 



without the need of any further calculation. Notice that 
for very small values of the dispersion, the accuracy of 
our calculations is dominated by precision errors, and 
these expressions for the error bounds do not apply. 



V. CONCLUSIONS 

In this paper, we develop a new method to construct 
basis sets from scar functions which are able to calcu- 
late very accurately the eigenstates of classically chaotic 
systems and keeping the size of the problem at very mod- 
erate sizes. The main idea is that introducing the rele- 
vant dynamical information into the basis elements has 
a twofold effect. First, the efficiency of the basis set is 
fostered, and second, its nature allows us an straightfor- 
ward identification of the underlying classical structures 
contributing to each individual eigenstate, thus provid- 
ing a good description of the quantum dynamics in a 
semiclassical sense. 

As an illustration we applied the method to a clas- 
sically chaotic quartic two-dimensional oscillator [see 
Eq. 0], that has been used as a benchmark in the field 
of quantum chaos. The performance of our method for 
this system was superb, since we have demonstrated that 
the method can be advantageously applied to the com- 
putation of excited states in a small energy window using 
very modest basis sets, whose sizes can be estimated from 
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the average participation ratio. Moreover, using a basis 
set of ~2500 elements, we have been able to compute the 
first low-lying ~2400 eigenstates with an accuracy, both 
in energies and wave functions, similar to that of other 
standard methods used in the past to study the same 
oscillator. 

Furthermore, we have examined the quality of our re- 
sults using a variety of different indicators, such as pro- 
jection on local (scar) representations, scar intensities, or 
participation ratios. The first two allow a direct descrip- 
tion of the different eigenstates in terms of semiclassical 
ideas, while the latter provides an idea of how our scar 
basis is much smaller than other conventional one. Fur- 
thermore, we also provided in this work upper bounds to 
the errors that can be expected in our calculations, both 
in energy and wave functions. 

We are currently extending the method to realistic 



molecular systems with a mixed phase space, where the 
degree of chaoticity depends on the energy. The results 
of an application to the LiNC/LiCN isomerizing system 
will be reported elsewhere. 
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SUPPLEMENTAL MATERIAL 

Here, we present some additional material containing 
our results on the structure of the eigenfunctions of the 



quartic oscillator in a more systematic way. The rele- 
vant numerical data corresponding to the reconstruction 
with scar functions of the eigenfunctions 1 220 + 5i)^ 15 
with i = 0,1,..., 23, whose probability densities are 
shown in Fig. [20j are reported in Table |VI[ In it, we 
have included all scar functions necessary to reconstruct, 
in each case, more than 85% of the corresponding eigen- 
functions, remark that all computed states were obtained 
with the same accuracy as obtained by other standard 
methods |TJ|2]. As we can see in the figure, the pattern 
of several eigenfunctions is concentrated only over a sin- 
gle PO of the system. In these cases, the first element of 
the local representation coincides with the scar function 
that is localized over that specific PO, similarly to what 
happens with the eigenfunctions in Figs. 9-11. In the 
rest of the cases, and although the eigenfunctions are in 
general not scarred by a single PO, the SGSM does still 
permit their reconstruction using a few (localized) scar 
functions of our basis set, as can be checked in the data 



presented in Table VI 



In a similar way, we present in Tables |VIIfjX| the re- 
sults corresponding to the eigenfunctions in Figs. 21 22} 
associated to the other symmetry classes. 
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Table VI. Reconstruction of the eigenstates \N)a ± with the scar functions |PO, ti)ai- The participation ratio, R, defined in 
Eq. (35) and the accumulated reconstruction overlap, £%, are also given. 

N R PO,n,E% 



220 2.37 
225 3.83 
230 2.71 
235 6.45 
240 3.45 
245 6.43 
250 3.83 
255 6.52 
260 3.56 
265 1.30 
270 4.60 
275 5.25 
280 2.15 

285 4.27 
290 6.87 

295 3.45 
300 2.50 
305 8.03 
310 2.44 
315 8.96 

320 3.60 
325 2.93 
330 9.38 

335 6.55 



14,44,62.5 
4,58,35.3 
14,45,57.6 
3,19,20.0 
14,46,49.8 
12,30,29.5 
16,40,46.7 
1,52,26.1 
9,52,49.7 
1,53,87.7 
9,53,36.6 
2,54,31.0 
9,54,67.7 
2,55,86.4 
2,55,41.0 
9,55,29.3 
13,36,74.7 
4,67,48.9 
9,56,58.9 
4,68,13.4 
9,57,62.8 
11,55,23.0 
16,45,74.2 
9,58,49.0 
16,46,52.1 
10,58,18.9 
3,23,84.5 
18,50,35.1 
13,39,84.6 



4,57,77.6 2,48,82.9 9,47,86.5 
9,48,63.0 6,28,86.7 

2,49,74.0 7,21,77.7 6,28,81.3 16,38,85.6 

4.59.43.5 17,36,61.3 7,21,76.3 14,46,79.1 14,45,81.3 16,39,82.4 13,33,83.4 5,12,84.4 6,28,85.1 

5.12.62.6 16,39,75.9 13,33,82.4 4,59,86.1 

4,61,42.5 17,37,52.9 6,29,58.4 16,40,63.4 1,51,66.8 2,51,82.0 9,50,86.1 

12,30,63.3 9,51,70.5 17,37,75.2 7,22,78.8 2,51,82.3 1,51,85.6 

4,62,43.9 14,48,61.0 3,20,72.3 13,34,78.9 7,22,85.7 

14,48,58.3 15,46,68.4 6,30,76.6 4,63,82.8 12,31,88.1 

4,64,59.2 17,39,72.1 8,31,80.5 13,35,87.4 

5,13,55.9 8,31,69.0 13,35,76.0 4,64,84.9 17,39,87.5 

11,52,70.9 17,40,72.9 13,36,75.5 10,53,76.4 5,13,77.2 10,54,77.9 9,55,79.4 4,66,80.9 1,55,81.9 



13,36,59.9 
4,66,39.6 
16,43,75.9 
10,55,60.9 
16,44,73.9 
18,48,26.9 
11,54,69.5 
2,58,41.3 
6,33,75.5 
11,55,57.3 
8,34,77.4 
4,71,33.0 
16,47,87.1 
17,44,42.8 
2,59,87.4 



3,21,71.6 17,40,78.6 4,65,88.7 

10.54.47.0 15,49,49.5 6,32,51.3 8,32,67.2 7,24,69.7 17,41,71.5 4,67,72.9 11,52,73.8 
1,56,76.6 2,56,86.7 

9,56,75.6 17,41,85.1 
10,55,90.8 

10,56,42.3 11,54,56.6 9,57,70.4 3,22,85.6 

6.33.73.2 18,48,76.0 10,56,77.4 1,58,78.7 2,58,86.1 

12.34.49.6 4,69,54.3 5,14,58.1 1,58,60.3 9,57,67.8 11,54,69.7 7,25,71.2 18,48,72.6 

3.22.76.3 10,57,77.0 9,58,80.3 18,49,86.9 

4,70,63.2 5,14,65.1 1,59,66.7 2,59,81.9 11,56,84.1 17,43,87.0 
17,43,82.8 7,25,86.3 

18.50.39.7 11,56,44.2 6,34,48.2 16,46,56.5 8,34,61.3 1,60,64.5 2,60,81.2 13,39,82.3 

11.57.50.1 4,71,54.4 2,60,58.4 10,58,62.8 1,60,67.0 11,56,70.7 16,46,74.4 6,34,80.1 



Table VII. Same as Table VI for the eigenfunctions \N) 



/A 2 - 



N R PO,n,£% 



220 4.67 6,28,40.5 4,58,49.7 7,21,62.9 9,47,68.5 2,48,79.7 11,45,88.1 

230 3.17 4,59,43.3 18,42,77.9 10,48,79.4 10,49,81.0 11,46,82.4 9,48,90.6 

240 4.62 14,47,43.8 10,50,53.1 4,60,60.2 5,12,64.0 6,29,66.8 17,38,68.9 7,22,72.9 12,31,76.2 16,40,79.8 2,51,83.3 
8,30,85.1 

250 4.76 8,30,37.7 4,62,56.9 17,39,65.3 10,51,73.8 6,30,84.4 16,41,88.8 

260 4.19 14,49,39.7 2,53,53.3 10,52,76.3 12,32,82.2 9,52,85.8 

270 4.13 16,42,42.7 13,35,49.8 17,40,59.9 14,50,79.4 4,64,85.2 

280 3.79 2,55,38.5 12,33,68.9 14,51,81.5 4,66,85.2 

290 1.49 2,56,81.5 4,67,86.1 

300 3.56 14,53,47.7 10,56,65.1 8,33,77.8 18,49,82.1 9,56,85.1 

310 1.75 2,58,72.2 13,38,94.4 

320 4.90 6,34,33.0 18,50,62.0 11,55,66.7 2,59,70.9 15,52,74.7 9,57,77.5 5,14,81.6 2,58,84.6 4,70,87.2 

330 1.22 2,60,90.5 




Figure 20. (Color online) Quartic oscillator eigenfunctions \N)a 1 corresponding to the reconstruction data of Table VI Same 
coordinates and scaling as in Fig. 5. 



24 



Table VIII. Same as Table VI for eigenfunctions \N)b 1 - 



N R PO,n,£% 



220 3.03 
230 5.20 
240 7.76 

250 2.36 
260 1.76 
270 2.35 
280 6.48 
290 4.40 
300 2.25 
310 8.11 

320 7.56 
330 7.09 



14,45,46.8 
4,59,31.4 
6,29,29.2 
13,33,76.3 
16,40,63.9 
9,52,73.5 
8,31,64.2 
5,13,25.9 
6,32,43.0 
13,37,63.3 
10,57,14.0 
6,33,80.5 
12,35,28.3 
10,59,25.5 
4,70,83.5 



2,49,79.4 
16,38,58.6 

9,50,45.3 
17,38,77.9 
11,49,71.3 
17,39,89.7 
18,46,70.1 

1,56,36.8 
13,36,58.3 

4,68,83.7 

1,59,19.1 
12,35,83.5 

4,70,42.7 

4,71,35.7 
17,43,84.3 



16,37,84.3 
17,37,66.7 
16,39,54.2 
4,61,79.3 
7,22,77.9 

12,32,76.7 
2,56,54.5 
4,67,65.9 
18,48,87.6 
2,59,47.7 
2,58,86.5 
13,38,53.9 
17,44,41.2 
11,57,85.0 



17,36,87.7 

7.21.73.0 6,28,76.8 2,50,79.1 1,50,82.2 9,48,89.3 

5,12,57.6 11,48,62.1 2,51,65.9 14,47,68.4 12,30,71.3 10,51,72.8 9,51,74.6 

7.22.81.8 4,60,83.2 2,52,84.6 12,31,85.8 
13,33,79.0 17,38,83.9 5,12,86.7 

4,64,81.2 9,53,83.7 14,51,85.2 

9.54.72.1 12,33,78.7 18,47,87.3 

5.13.72.9 11,52,74.5 4,66,76.2 11,53,77.4 1,57,78.6 2,57,83.8 9,55,88.4 

9,57,54.8 11,55,61.2 11,54,65.0 4,69,68.3 17,43,72.2 13,37,74.7 18,49,77.2 

18,50,59.9 8,34,65.4 5,14,71.4 6,34,75.8 1,60,78.0 2,60,81.9 9,58,86.8 

1,61,46.0 2,61,69.4 5,14,73.8 8,34,75.5 6,34,79.4 12,35,81.1 13,38,82.2 



Table IX. Same as Table VI for eigenfunctions \N)b 2 



N R PO,n,S% 



220 7.55 10,46,27.1 4,57,32.8 9,47,36.9 2,47,51.7 11,45,67.2 13,32,70.0 10,47,72.1 6,28,73.5 8,28,79.0 3,19,80.8 

2,48,82.1 5,12,83.4 11,46,84.1 9,48,88.9 

230 4.30 18,41,37.4 5,12,64.4 10,47,75.7 8,28,76.8 9,48,77.8 11,46,84.1 2,48,87.6 

240 1.53 4,60,79.7 18,42,92.8 

250 2.72 18,43,50.0 14,47,83.6 16,40,88.5 

260 2.00 14,48,70.3 11,50,73.5 9,51,75.0 5,13,76.4 10,50,77.5 2,51,80.5 10,51,82.3 14,47,83.4 4,62,84.4 11,49,86.5 

270 3.26 4,64,47.4 6,31,71.7 14,49,82.5 16,42,92.9 

280 3.93 16,43,44.0 17,40,61.8 10,53,76.1 2,54,82.9 14,50,84.9 9,54,89.2 

290 3.79 14,51,46.7 17,41,60.7 4,67,67.6 15,49,75.0 7,24,86.2 

300 8.66 14,52,22.3 2,56,39.4 6,33,47.1 10,55,53.5 9,56,62.8 4,68,68.9 3,22,76.1 18,48,78.8 15,50,83.4 5,14,90.3 

310 3.71 12,34,45.0 4,69,67.3 7,25,76.5 15,51,83.0 6,33,87.6 

320 4.58 13,39,32.9 6,34,61.3 10,57,69.2 2,58,74.8 12,35,79.9 14,54,82.9 9,58,95.7 

330 4.02 16,47,46.2 3,23,55.2 17,44,63.6 10,58,72.2 7,26,77.8 18,50,84.4 13,40,90.5 




Figure 21. (Color online) Same as Fig. 20 for the eigenfunctions \N)a 2 and \N)b 1 corresponding to the reconstruction data of 
Tables IVlTI and IVTTTl 
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Table X. Same as Table VI for eigenstates \N)e 1 - 



N R PO,n,£% 



220 4.47 
230 4.53 

240 4.34 
250 2.15 
260 2.85 
270 4.35 
280 6.75 

290 6.35 
300 6.94 
310 4.11 

320 6.70 
330 9.24 



9,34,40.2 

10,32,42.9 
9,35,84.8 

15,30,30.3 
2,37,62.3 
4,43,52.8 

12,22,40.3 
4,45,31.4 
9,38,83.8 

15,35,28.3 
9,40,27.2 

15,36,47.3 
4,48,82.1 
4,50,28.0 
4,49,19.3 

11,40,85.0 



7,27,54.3 
13,23,55.6 
9,34,85.8 
2,36,58.4 
15,31,89.5 
12,22,78.3 
4,44,62.4 
7,31,46.4 
15,34,84.7 
6,45,43.1 
7,32,45.4 
12,24,55.3 
7,33,83.1 
3,16,41.1 
12,25,29.6 



2,31,67.6 14,32,81.8 8,39,86.2 

12.20.66.0 2,35,73.7 4,40,77.4 8,39,79.0 15,29,79.8 6,40,80.7 4,42,82.4 11,33,83.8 

11,34,81.7 5,17,85.7 

5.18.83.6 4,44,87.4 

6,43,70.8 13,25,76.3 10,35,81.8 8,44,86.8 

2,36,56.1 12,23,64.4 14,34,69.9 10,36,76.8 9,39,80.1 10,38,31.4 3,15,82.5 2,40,83.3 
11,38,85.2 

4,46,61.4 12,23,73.2 3,15,76.2 13,26,83.2 11,38,85.9 

12,23,60.2 4,46,66.0 15,35,69.3 2,37,73.0 4,48,75.6 11,38,77.8 10,37,85.0 

11,39,60.2 6,46,63.0 8,46,71.3 10,40,43.5 9,40,75.2 14,36,76.9 9,38,79.0 12,23,80.8 

15.35.84.1 7,32,84.5 13,27,84.8 4,49,85.2 
10,41,55.4 15,37,66.1 2,39,76.9 2,42,85.7 

4.50.37.7 6,20,46.2 13,28,63.1 3,16,75.1 4,51,78.3 9,40,79.7 2,43,82.0 14,37,83.2 
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